https://github.com/satellogic/telluric-talks




Resumen de proyectos clave:
| Python | C/C++ | |
|---|---|---|
| pyproj | ⇒ | PROJ.4 |
| Fiona | ⇒ | OGR |
| Shapely | ⇒ | GEOS |
| rasterio | ⇒ | GDAL |
Ventajas:
Desventajas:
telluric is an open source (MIT) Python library to manage vector and raster geospatial data in an interactive and easy way.
Importar para uso interactivo:
import telluric as tl
from telluric.constants import WGS84_CRS, WEB_MERCATOR_CRS
Basic geometry definition using GeoVector: in the simplest case, the bounds and the projection (CRS)
ebro = tl.GeoVector.from_bounds(
xmin=0, ymin=40, xmax=1, ymax=41, crs=WGS84_CRS
)
print(ebro)
GeoVector(shape=POLYGON ((0 40, 0 41, 1 41, 1 40, 0 40)), crs=CRS({'init': 'epsg:4326'}))
ebro
También se pueden usar geometrías de Shapely:
from shapely.geometry import Polygon
ebro2 = tl.GeoVector(
Polygon([(0, 40), (1, 40.1), (1, 41), (-0.5, 40.5), (0, 40)]),
WGS84_CRS
)
ebro2
ebro.centroid
ebro.area # Real area in square meters
9422706289.175217
ebro.within(ebro2)
False
print(ebro.difference(ebro2))
GeoVector(shape=MULTIPOLYGON (((0 40.66666666666666, 0 41, 1 41, 0 40.66666666666666)), ((1 40.1, 1 40, 0 40, 1 40.1))), crs=CRS({'init': 'epsg:4326'}))
ebro.difference(ebro2)
GeoVector + atributosFeatureCollections: una secuencia de featuresgf1 = tl.GeoFeature(
ebro,
{'name': 'Delta del Ebro'}
)
gf2 = tl.GeoFeature(
ebro2,
{'name': 'Delta del Ebro + Castellón'}
)
print(gf1)
print(gf2)
GeoFeature(Polygon, {'name': 'Delta del Ebro'})
GeoFeature(Polygon, {'name': 'Delta del Ebro + Castellón'})
fc = tl.FeatureCollection([gf1, gf2])
fc
# This will only save the URL in memory
rs = tl.GeoRaster2.open(
"https://github.com/mapbox/rasterio/raw/master/tests/data/rgb_deflate.tif"
)
# These calls will fetch some GeoTIFF metadata
# without reading the whole image
print(rs.crs)
print(rs.band_names)
CRS({'init': 'epsg:32618'})
[0, 1, 2]
rs.footprint() # Bounding box
rs
rs.image[:,200:300, 200:240]
masked_array(
data=[[[15, 20, 23, ..., 41, 38, 29],
[17, 15, 23, ..., 41, 46, 46],
[18, 17, 17, ..., 39, 39, 44],
...,
[106, 49, 62, ..., 46, 58, 50],
[255, 255, 255, ..., 38, 37, 43],
[255, 229, 255, ..., 40, 31, 35]],
[[94, 98, 100, ..., 129, 98, 62],
[98, 96, 100, ..., 127, 127, 133],
[98, 96, 94, ..., 126, 126, 129],
...,
[145, 79, 116, ..., 47, 58, 51],
[255, 255, 255, ..., 43, 41, 45],
[255, 249, 255, ..., 41, 33, 37]],
[[131, 131, 131, ..., 161, 125, 86],
[133, 133, 135, ..., 164, 160, 161],
[133, 132, 129, ..., 159, 160, 157],
...,
[141, 77, 108, ..., 29, 34, 31],
[255, 255, 255, ..., 31, 30, 31],
[255, 255, 255, ..., 23, 19, 27]]],
mask=[[[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
...,
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]],
[[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
...,
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]],
[[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
...,
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]]],
fill_value=999999,
dtype=uint8)
# Crop fetch only the data required to build the raster
rs.crop(rs.footprint().buffer(-50000))
rs[200:300, 200:240] # the image is streched due to the presentation layer
# Rasterizing Vectors
# Geting OpenStreetMap data of Tel Aviv
fc = tl.FileCollection.open("data/telaviv.geojson")
fc
streets_raster = fc.rasterize(dest_resolution=1)
streets_raster
.crop y .rasterization para distribuir la carga en varias máquinasfrom telluric.vectors import generate_tile_coordinates
list_of_regions = list(generate_tile_coordinates(rs.footprint(), num_tiles=(10,10)))
# This line is to visualize and not required
tl.FeatureCollection.from_geovectors(list_of_regions)
chunk0 = rs.crop(list_of_regions[17])
chunk0
import dask.multiprocessing
from dask.delayed import delayed
from dask.base import compute
def worker(raster_filename, roi):
raster = tl.GeoRaster2.open(raster_filename) #Lazy loading of the raster
chunk0 = raster.crop(roi)
return chunk0.image.max()
RASTER_URL="https://github.com/mapbox/rasterio/raw/master/tests/data/rgb_deflate.tif"
items = [delayed(worker)(RASTER_URL, roi) for roi in list_of_regions]
maxs = compute(*items, get=dask.multiprocessing.get)